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I.  INTRODUCTION 


This  report  is  concerned  with  least-squares  tracking  algorithms  predicted  on  mini¬ 
mizing  frequency  and  bearing  errors.  The  state  variables  which  we  wish  to  determine  may  be 
considered  to  be  the  initial  target  position  (x-coordinate  and  y-coordinate)  and  target  velocity 
for  each  segment  of  the  target  track  (2(  number  of  maneuvers  +  1 )  variables).  The  bearing 
angles  0r  and  frequencies  f*  will  be  a  function  of  the  unknown  state  variables  q  =  (qj.  q2- 
■  •  •)  as  well  as  of  parameters  associated  with  the  receiving  platforms  (position,  velocity, 
etc.)  which  we  denote  by  the  vector  z.  The  superscript  “t”  is  intended  to  suggest  time,  but 
can  actually  index  any  general  set  of  measurements  not  necessarily  in  chronological  order 
and  possibly  received  by  several  platforms. 

In  terms  of  the  above  definitions,  the  problem  may  be  expressed  as: 

Determine  q  which  minimizes 


T  T 

G(q)  =  ^  at(0l 2 3 4(q,  z)  -  ^ 

t=l  t=l 


fl  (q,  z) 


(1) 


where  0**  and  T*  are  measured  quantities  and  a{  and  bt  are  suitable  weights  (set  equal  to 
zero  in  the  absence  of  a  measurement).  Note  that  the  source  frequency  fQ  can  either  be 
considered  a  parameter  (such  as  z)  or  one  of  the  state  variables  q.  A  necessary  condition 
for  q  to  yield  a  minimum  is  that  it  satisfy  the  set  of  nonlinear  equations 


VqG(q) 


a  =  0 
q 


(2) 


obtained  by  setting  the  gradient  of  ( 1 )  equal  to  zero. 

An  obvious  choice  for  the  weights  is  to  make  them  proportional  to  the  accuracy  of 
the  corresponding  measurements.  If  we  take  at  =  and  bt  =  (o£t  j-1  where  a~ 


denotes  variance  we  obtain  the  so-called  maximum  likelihood  estimator  ( MLE)  [1,2].  This 
nomenclature  reflects  the  fact  that  if  (T1  and  are  independent  Gaussian  random  variables. 

A.  _ 

then  the  minimum  q  is  the  maximum  likelihood  estimate  of  q.  Bowing  to  current  usage 
[2] ,  we  shall  continue  to  use  the  term  MLE  although  it  is  a  misnomer  when  at  and  bt  are 
not  the  true  variances. 


There  are  many  standard  techniques  for  solving  (2).  We  merely  note  here  that  they 
typically  involve  the  computations  of  the  gradient  and/or  the  Jacobian  of  G([2] ,  [3] ,  [4] ). 
These  computations  are  discussed  in  Appendix  A.  Our  primary  concern  in  this  report  is 


1.  Jazwinski,  A.,  Stochastic  Processes  and  Filtering  Theory.  Academic  Press,  N.Y.,  1970. 

2.  Adams,  W.B.,  “A  Non-Linear  Maximum  a  posteriori  Estimator  for  Passive  Ranging  Using  Bearing  and 
Frequency  Measurements  on  a  Stable  Line,”  General  Electric  Co.,  Syracuse,  N.Y.,  Sep  1975. 

3.  Powell,  M.,  “A  Survey  of  Numerical  Methods  for  Unconstrained  Optimization,”  SIAM  Review.  Vol  12. 
pp  79-97,  Jan  1970. 

4.  Bard,  Y„  “Comparison  of  Gradient  Methods  for  the  Solution  of  Nonlinear  Parameter  Estimation 
Problems,”  SIAM  J.  Num.  Analysis,  Vol  7,  pp  157-185,  March  1970. 
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the  accuracy  of  the  solution  if  as  a  function  of  the  statistics  of  <F.  f  .  and  z.  This  is 
treated  in  the  next  section.  The  analysis  permits  the  direct  calculation  of  the  error  statistics 
(to  a  reasonable  approximation),  a  capability  which  is  important  in  practical  implementa¬ 
tions  and,  as  well,  alleviates  the  need  for  Monte  Carlo  techniques  in  evaluations  of  the 
algorithm. 

We  shall  also  see  that  for  Gaussian  statistics,  the  estimated  MLE  error  is  essentially  the 
Cramer-Rao  bound  [5J  for  the  solution  track.  If  the  errors  are  not  too  large  this  estimate 
a  good  approximation  of  the  actual  error  statistics:  i.e.,  the  MLE  achieves  minimum  variance 
compared  to  all  unbiased  estimators  (with  white  Gaussian  noise).  A  similar  result  has  been 
shown  for  the  Kalman  filter*  linearized  about  the  true  solution  16) .  Thus,  if  one  assumes 
the  bias  is  negligible  (as  it  appears  in  practice),  the  local  accuracy  of  the  MLE  and  extended 
Kalman  filter  should  be  approximately  the  same,  and  the  error  analysis  developed  here 
should  serve  as  a  lower  bound  for  this  entire  class  of  algorithms. 

Furthermore,  the  computations  in  the  following  section  include  a  sensitivity  analysis 
with  respect  to  most  of  the  significant  measured  parameters  such  as  platform  position,  speed, 
etc.  (designated  z  in  equation  ( 1 )).  The  validity  of  these  results  extends  to  any  algorithm 
designed  to  minimize  the  MLE  cost  function,  and,  thus,  should  also  be  indicative  for  the 
various  nonlinear  Kalman  filters. 

In  brief,  we  shall  first  derive  general  expressions  for  the  error  by  linearizing  and  taking 
expectations.  These  are  then  specialized  to  cost  functions  of  the  type  ( 1 ).  Specifically,  we 
distinguish  between  the  roles  of  the  random  variables  z  and  those  of  6  and  f.  Finally,  we 
show  that  this  procedure  leads  to  an  estimate  of  the  variances  of  the  state  variables  equal  to 
the  Cramer-Rao  bound.  Appendices  A  and  B  contain  explicit  computations  of  equation 
(2).  the  Jacobian  of  that  system  of  equations  evaluated  at  the  solution  point,  and  the 
covariance  matrix  for  the  state  variables.  Appendices  B  and  C  also  present  a  fairly  detailed 
analysis  of  the  behavior  of  the  corresponding  confidence  ellipsoids.  Various  remarks  inter¬ 
preting  the  theory  and  pointing  out  its  (sometimes)  subtle  consequences  are  scattered 
throughout  the  text  and  are  considered  to  be  one  of  the  main  contributions  of  this  report. 

In  particular,  the  examples  found  in  Appendix  B  should  not  be  overlooked. 


*For  linear  systems,  the  Kalman  filter  (with  unknown  initial  state)  conicides  with  the  MLE  ( 1 1 :  however, 
nonlinear  versions  such  as  the  Iterated  Extended  Kalman  Filter  (1EKF)  are.  in  a  sense,  only  approximations 
to  the  MLE.  The  Kalman  filter,  being  a  recursive  algorithm,  generates  new  states  utilizing  only  the  previous 
state  and  the  most  recent  observation.  For  linear  estimators  this  is  sufficient  since  the  “linear  least  squares” 
or  “minimum  variance”  correction  is  a  linear  combination  of  the  observation  and  old  state  (orthogonal  to 
the  old  state).  For  nonlinear  estimators,  the  new  estimate  is  not  necessarily  orthogonal  to  the  old  since  the 
new  estimation  space  is  entirely  different:  it  involves  nonlinear  combinations  with  past  data  and  is.  in  genera), 
infinite  dimensional.  In  many  situations  the  approximation  involved  in  the  IEKF  are  good:  however,  they 
can  diverge,  and  the  IEKF  covariance  estimate  is  often  overly  optimistic  [1.7], 


1.  Jazwinski,  A.,  Stochoastic  Processes  and  Filtering  Theory,  Academic  Press.  N.Y..  1970. 

5.  Van  Trees.  H.L.,  Detection.  Estimation,  and  Modulation  Theory.  Part  I.  Wiley.  New  York.  N.Y.,  1968. 

6.  Taylor,  J.,  “The  Cramer-Rao  Estimation  Error  Lower  Bound  Computation  for  Deterministic  Nonlinear 
Systems,”  Proceedings  of  1978  IEEE  Conference  on  Decision  and  Control,  pp  1 178-1 181 . 

7.  Chou,  S.I.,  “Some  Drawbacks  of  Extended  Kalman  Fillers  in  ASW  Passive  Angle  Tracking,”  ONR 
Symposium  on  Advances  in  Acoustic  Passive  Tracking,  NPS-62  TS  77771 ,  Naval  Post  Graduate 
School,  Monterey,  CA.  May  1977. 
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2.  SENSITIVITY  AND  ERROR  ANALYSES 

Let  q  be  a  vector  representing  the  set  of  state  variables  (target  position,  velocity,  etc.) 
and  let  w  be  a  (scalar)  measured  parameter  such  as  receiver  position  or  target  bearing 
measurement.  Assume  that  we  are  solving  a  system  of  equations  in  q  of  the  form 


Fk(q,w)  =0  k  =  1,  •  •  -,N  • 

(3) 

The  sensitivity  to  the  measurement  w  may  be  computed  by  varying  w  from  its  true  value  w. 

A  A  A 

and  linearizing  about  the  point  (q,  w)  satisfying  Fk(q,  w)  =  0. 

That  is. 

Y  (Aq)m+  3[>Aw~0 

m  3qm  '  3w 

(4) 

or 

A’~-r'  (i) Aw  ■ 

(5) 

where 

A 

Aq  =q  -q 

(6a) 

A  w  =  w  -  w 

(6b) 

0Fk 

Jkm=  0^,  ,A  A 

m  (q,  w) 

(6c) 

and  F  is  an  N-dimensional  vector  valued  function,  with  components  Fk  . 

More  typically  we  are  concerned  with  functions  which  are  a  sum 
form  (3) 

of  terms  of  the 

^  F*.  (q,  wl)  =  0  k  =  l . N  . 

t 

(7) 

Let  us  assume  that  Aw^w-w*  for  t=l,  •  •  *,  T,  constitutes  a  set  of  zero-mean  independent 
identically  distributed  random  variables  such  as,  for  example,  the  set  of  bearing  measure¬ 
ments  at  time  t.  Then  we  may  write  E(AwlAwl  )  =  0  and  =  E(Awl)  .  Taking  the  expec¬ 

tation  of  the  product  of  components  of  (5),  we  obtain 

7 


(8a) 


T 


Let  us  now  assume  that  equation  ( 7)  is  the  result  of  a  nonlinear  least-squares  fit  (cf.. 
eq.  ( 1));  i.e.,  arises  from  the  minimization  problem 


min 

q 


(9) 


where  for  simplicity  of  notation  we  take  zx  to  be  a  scalar.  For  the  general  case  see  Appendix 
B  and  equation  (B-8).  Both  z£  and  <T£  represent  measured  parameters  ( for  example,  receiver 
speed  and  target  bearing  normalized  by  the  least  square  weights  respectively)  and  thus  are 
examples  of  the  variable  w£  in  (.7).  However,  their  roles  are  typically  perceived  as  being 
physically  different,  and  they  are  distinguished  mathematically  by  the  different  functional 
dependencies  of  expression  (9).  More  intuitively,  we  might  say  that  the  c"*  are  those  para¬ 
meters  which  have  been  chosen  to  "measure”  the  quality  of  fit  of  the  least  squares  track. 

The  solution  to  (9)  satisfies  the  set  of  equations 


I 


(c* 


c  l> 


dc^ 

3qk 


=  0 


k=l.  •  •  •,  N 


(10) 


Noting  that  at  the  point  of  linearization  (q,z£)  and  for  the  true 
the  expression  c£  -  c  £  vanishes,  we  find  the  Jacobian  of  system  ( 1 0)  to 


measured  value  c  , 
be  given  by 


,  .A  a  £ . 

Jkm(q.z.c  ) 


S  3c£  3c^ 
^  —  * 
t  3qk  3qm 


Also, 


3c* 

A  .  A 


(ID 


(12a) 


8 


.  =  -  p  b* 

lo 


We  note  that  the  lust  equation  is  equivalent  to 


=  -  —  .  All  derivatives  were 


evaluated,  of  course,  at  the  point  representing  the  true  parameter  values. 

Finally,  in  either  case,  if  we  are  interested  in  the  errors  due  to  the  measurements  f  l. 


°f  A  °fl 


tire  appropriate  relationship  is  ( 8)  with  o~,  replaced  by  —  =  — —  and  —7^—  given 

*•-  And 


(l  2b)  or  (16). 


REMARKS 

1.  Note  that  the  entire  analysis  rests  on  the  approximation  (4).  Although 
experience  would  lead  one  to  believe  that  (4)  is  a  valid  step  (  for  sufficiently  small  ax)  in  the 
derivation  of  (8).  one  must  still  retain  a  degree  of  skepticism.  For  example,  taking  the 
expectation  of  (4)  yields  F<Aqm)  =  0:  i.e..  seems  to  imply  that  the  least-squares  estimates 
are  unbiased.  This  is  only  true  to  the  first  order:  and.  as  is  well  known,  for  many  cost 
functions  the  bias  can  be  significant  < cf.  (8  j  and  the  Churn  algorithm  of  [7| ). 

2.  Under  the  assumption  that  the  various  types  of  measurement  errors  are 
independent  of  each  other,  we  may  calculate  E(Aq  )-  due  to  all  these  sources  simply  by 
summing  the  individual  contributions  (sec  Appendix  B).  The  RMS  error  is  then  found  by 
taking  the  square  root. 

3.  Relation  (ID  which  expresses  the  Jacobian  at  the  “true  solution"  point  holds 

at  the  solution  for  any  minimization  problem  of  the  form  min  £(G1)-  provided  the  minimum 
is  zero. 

4.  The  above  results  easily  generalize  to  cost  functions  which  are  the  sum  of 
functions  of  the  form  (4);  e.g.. 


mi"  \  Y  !  (c‘  -  c  >r  +  <d<  -  dV  1 
q  2 

t 


7.  Chou,  S.I.,  “Some  Drawbacks  of  Extended  Kalman  falters  in  ASW  Passive  Angle  Tracking,"  ONR 
Symposium  on  Advances  in  Acoustic  Passive  Tracking.  NPS-P2  TS  77771 .  Naval  Post  Graduate  School. 
Monterey,  CA.  May  1977. 

8.  Eykhoff,  P.,  System  Identification,  Chapter  b.  John  Wiley  &  Sons,  N.Y..  1974. 
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T1 


where  c  and  d  represent  bearing  and  Doppler  measurements  respectively.  Also,  nowhere  has 
our  analysis  restricted  us  to  measurements  made  from  a  single  receiving  platform.  However, 
if  the  measurements  have  different  variances,  additional  terms  of  the  above  form  might  be 
desirable. 


CRAMER-RAO  BOUND 

Let  us  assume  that  the  least  squares  weights  have  been  chosen  equal  to  the  inverse 
variances  of  the  measured  quantities  (true  MLE)  so  that  the  variances  of  the  measured 
parameters  of  expression  (9)  are  equal  to  one.*  If  these  variables  are  indepen¬ 

dent  Gaussian  variables,  the  likelihood  function,  up  to  a  constant,  is  given  by 

L(q)  =  -  ^  ^  -  c"1)  2  (17) 

t 

The  Fisher  information  matrix  is  (5) 

4,= 


Thus,  the  Cramer-Rao  bound,  which  is  true  for  any  unbiased  estimate  of  qk  (5] ,  may  be 
expressed  as 

E(Aqk)2  >  (r')kk  (19, 

On  the  other  hand,  our  sensitivity  analysis  led  to  equation  (8b),  which  with  o ^  = 

=  1,  gives 
c  1 


l  (s '  -  c~<)  ♦  y  *  * 


3c*  3cl 


3qk  3q 


m 


Jkm 


(18) 


*As  stated  previously,  these  parameters  are  also  assumed  independent.  The  more  general  case  may  be 
treated  by  using  the  inverse  covariance  matrix  in  lieu  of  the  weights. 

5.  Van  Trees,  H.L.,  Detection,  Estimation,  and  Modulation  Theory,  Part  I,  Wiley,  New  York,  N.Y.,  1968. 


II 


k 


r1  — 

km  aq 


t\  2 


m  i 


(from  1 2b) 


V  1  1  _dc^_  3c1 

^  Jkm  Jkm  3qm  9qm. 

m' 


-  2  2, 

m  m' 


j- 1 

Jkm 


rkln'  ' 


mm 


—  [—  1 
"  Jkk 


(20) 


Thus,  in  this  case,  our  estimated  error  of  the  MLE  algorithm  equals  the  Cramer-Rao  bound. 
The  same  computation  applied  to  (8a)  also  yields 

Q  =  J-1  (21) 


where  0  is  the  covariance  matrix  of  Aq. 

Note  that  each  new  data  point  results  in  the  addition  of  a  matrix  of  the  form  C]^ 
0ct  0^,t 

=  t=—  ^=7—  .  Since  C  is  positive  semi-definite,  J  is  a  non-decreasing  matrix  function 

3qk  3qm 

(J  >  J'  if  J— J '  is  positive  semi-definite)  of  the  data.  We  conclude  from  Appendix  C  that  the 
corresponding  concentration  ellipsoids  (Aq)  ^Q(Aq)  <  1  (Appendix  B  and  [51  )form  a  non¬ 
increasing  nested  sequence.  In  practice  C  is  typically  positive  and  the  sequence  is  strictly  de¬ 
creasing.  The  reader  should  be  warned,  however,  that  although  this  represents  an  improved 
estimate  of  the  state  variables,  it  does  not  necessarily  imply  increased  accuracy  in  the  current 
target  position.  The  position  state  variables  are  the  coordinates  of  the  target  at  a  fixed  point 
in  time  (for  example  t=0),  and  the  errors  must  be  propagated  along  the  target  track  as  de¬ 
tailed  in  Appendix  B. 


5.  Van  Trees,  H.L.,  Detection,  Estimation,  ami  Modulation  Theory,  Part  I,  Wiley.  New  York.  N.Y.,  1968. 


A  second  comment  is  that  relations  (20)  and  (21)  assume  that  the  a  priori  estimates 
of  the  measurement  variances  are  correct;  i.e.,  at  =  1  ./o$  and  bt  =  f q/o^  in  equation  ( 1 ). 
If  this  is  not  true, 


*K)2  -in  ^ m  2  a  afe 


t  m  m 


m 


and  further  reduction  is  not  generally  possible.  These  comments  are  illustrated  in  figures 
la  and  lb.  Inclusion  of  the  errors  due  to  the  parameters  z*  of  (9)  could  have  similar  effects; 
in  particular,  Q  #  J"  * . 


Figure  la.  Example  illustrating  diminishing  error  ellipses.  Platform  track  course  is 
20  minutes  at  0°  and  1 2  knots;  and  20  minutes  at  15°.  10  knots.  Target  is  0°.  25 
knots.  Bearing  standard  deviation  is  Og=0.4°  and  frequency  standard  deviation  is 
CTf=0.01  Hz  with  a  base  frequency  of  100  Hz.  The  MLE  weights  were  taken  to  be 

a~Q~  and  (of/f0)"“.  respectively.  (The  first  two  points  have  infinite  error  ellipses.) 


_2 

Figure  lb.  The  same  as  Figure  la.  only  the  MLE  weights  were  taken  to  be  Oq  and 
(0.5/fQ)  “  which  corresponds  to  a  Of  of  0.5  Hz.  Note  that  only  the  ratio  of  the  two 
weights  is  relevant. 


14 


APPENDIX  A 

COMPUTATION  OF  GRADIENT  AND  “JACOBIAN” 


Letp(t)  =  (px(t),  Py(t))  and  v(t)  =  (vx(t)  vy(t))  represent  the  position  and  velocity 

of  one  of  the  receiving  platforms  at  time  L  Similarly  let  p'(t)  and  v'(t)  be  those  of  the  source 
(target).  Their  relative  position  and  velocity  are  given  by: 

Ap  =  p'  -  p 

Av  =  v  -  v  (A-l) 

so  that,  if  we  define  the  positive  y  axis  as  north,  the  relative  bearing  of  the  target  is  given  by 
(cf.  Figure  A.  1) 


6  =  tan-* 


(Ap)x 

<Ap)y 


( A-2) 


(a) 


lb) 


Figure  A.l.  Example  of  the  geometry  for  a  two  segment  target. 


It  is  convenient  (and  provides  improved  roundoff  characteristics)  to  deal  with  the 
normalized  Doppler  shift 

f  -  fQ 

d  =  — -  -  (A-3) 

ro 


where  fQ  and  f  are  the  source  and  received  frequencies  respectively.  If  the  velocities  are 
small  with  respect  to  that  of  sound  (denoted  c),  we  may  consider  only  the  relative  motion 
in  computing  d.  This  results  in  the  following  expression: 
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d  =  - 


Av  •  Ap 

«ll  ®ll 


(A-4) 


where  “  *  ”  indicates  scalar  product. 

We  now  assume  that  the  target  follows  a  piecewise  linear  track  with  initial  position 
p^  and  velocity  vj  over  segment  i  from  time  tj  to  tj+j .  Then 

p(t)  =  %  +  vj  tt2~  t,>  +  +  ^(t-ti)  (A’5) 

for  tj^  t<ti+1 

.  .  v  , - '  v  <■£,'  v  ( (v'\  where  i=  l  to  the  number  of 

Our  state  variables  are  taken  to  be  ip0*X'  ^07’  vi  X’ '  1 7 

segments.  The  base  frequency  1'0  may  also  be  a  state  variable. 

The  cost  function  to  be  minimized  takes  the  form 


G 


JL  B  +  4  D 


°S 


°a 


( A-6) 


where 


(A-7) 


D  =  T  \ 


( A-8 ) 


and  a\  and  are  the  a  priori  variances  <  f_f  f- 

).  The  symbol  5”  represents  measured  bearing  and  3"  fQ 


riances  of  0  and  d  (for  notational  simplicity  assumed  con- 

fo 

_  .  .  _ _ urine  and  d  =  7“ 

stant  functions  of  t). 

whKe  f  is  ,he  measured  frequency.  The  base  frequency  fc  is  either  detemmed  a  priori 
f  nm  fhp  data  or  is  a  state  variable  to  be  included  in  the  minimization.  The  -ir  in  equation 
[a-7,  ind,"  aKs  fha,  the  bearing  difference  is  taken  modulo  >  so  as  to  yield  the  mm, mum 
difference,  i.e.,  a  value  between  -it  and  n. 

The  partial  derivatives  of  B  and  D  with  respect  to  a  state  variable  q  are 


(B)q  =  l  (*'-?')2*('X)  q 


lb 


T 


>D>q  ‘I 

t 


(A-9) 


with  the  subscript  q  indicating  partial  differentiation.  These  expressions  may  be  evaluated 
using  the  chain  rule.  Thus,  dropping  the  superscript  t  for  simplicity,  we  have 


(d  -  d)q 


|  (0)_,  •  (p')q  q  *  fG 

|  0  q  =  f0 


(d)_,  *  (p)fl  +  (d)  •  (Av)  q^fQ 

p  q  Av  q  u 


q=  fn 


( A- 1 0) 


( A- 1 1 ) 


The  stands  for  scalar  product;  a  vector  subscript  indicates  the  gradient;  and  the  deriva¬ 
tive  of  a  vector  is  the  vector  comprised  of  the  derivatives  of  its  components. 

The  gradient  vector  (0)^'  is  easily  found  from  (A-2). 


(0), 


(0), 


(Ap)y 

II  Apll 2 

-(Ap)x 

IIApIl2 


Likewise,  ( A-4)  implies 
-1 


(d)s'  = 


Px  IIApIl 
(Ap)x 


L 


(Ap)x  (Av), 


IIApIl 


(dV  = 


x  c||Apll 


(A- 12) 


C  A- 13) 


and  similarly  for  the  y  components. 

The  other  partials  follow  immediately  from  (A-5).  Some  of  these  values  depend  on 
the  time  t.  The  only  “tricky”  calculation  is 
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0 


3Px 

*ty* 


J  t  -  ti<t<tj+1 


(A- 14) 


We  also  note  that  (Av)q  is  zero  unless  q  =  (v-)x  or  (vj)y  where  tj  <  t  <  tj+| ,  in  which  case  it 
is  the  vector  ( 1 ,0)  or  ( 0, 1 )  respectively. 

The  computation  of  J  of  equation  ( 1 1 )  is  given  by 


( A- 15) 


Note  that  this  is  the  true  Jacobian  of  the  system  of  equations VG  =  0  only  at  a  point  for 

which  0l  =  6  *  and  dl  =  J1:  i.e.,  only  at  a  true  solution  with  perfect  measurements.  Since, 
in  practical  situations,  this  is  never  the  case,  J  only  represents  a  (possibly  poor)  approxima¬ 
tion  to  the  Jacobian.  Nevertheless,  it  is  quite  suitable  for  use  in  the  steepest  descent  type 
algorithms  employed  in  solving  (2)  [2). 

Finally,  we  remark  that  if  the  various  measurements  have  different  accuracies,  i.e., 
if  var(0  t)  is  a  function  of  t,  then  Oq  and  must  be  replaced  by  ofiit)  and  oy  t),  which 
reside  inside  the  summation  over  t. 


2.  Adams,  W.B.,  “A  Non-Linear  Maximum  a  posteriori  Estimator  for  Passive  Ranging  Using  Bearing  and 
Frequency  Measurements  on  a  Stable  Line,”  General  Electric  Co.,  Syracuse,  N.Y.,  Sep  1975. 
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APPENDIX  B 

STATE  VARIABLE  COVARIANCE  AND  ERROR  ELLIPSES 

The  equations  developed  in  Section  2  suffice  to  compute  the  covariance  matrix  of 
the  state  variables  p^  and  v-;  however,  we  are  typically  interested  in  the  accuracy  of  the 
estimated  current  target  position  p'(t)  (see  equation  (A-5)).  The  covariance  of  p'(t)  follows 
directly  from  (A-5) 


E(Apr(t)  Apq(t))  =  E  (Apor  +  ^T  Avir  Atj)  (Apoq  +  ^  Aviq  Atj) 


(B-l ) 


where  r  and  q  stand  for  either  x  or  y,  and  Atj=t-tj  or  tj  +  j  -  tj.  Let  the  covariance 
matrix  of  the  state  variables  be  given  by 

=  E(PorPoq) 


rq 


E(porviq) 


rq 


^  E(VirVjq) 


(B-2) 


Then  (B-l)  becomes 


E(Ap^tl  Apq(.»  =  Off  ♦  X  Ali  +  S  Atj  ^rqVj 


J 


(B-3) 


Note  that  the  matrices Q  may  be  computed  from  (8a),  (see  also  equation  (D-8)). 
Let  us  label  the  three  quantities  evaluated  in  (B-3)  by 

ffxx  =  E(Apx  Apx) 


Oyy  =  E(APy  A^y) 


Oxy  =  E(APX  APy)  • 

The  equation  of  the  ellipse  for  the  position  variables  is  then  (10] 
°-f  (Ax)*  -  (Ax)  (Ay)  +  Ai  (Ay, 2  =  c 


(B-4) 


(B-5) 


10.  Morrison,  D.,  Multivariate  Statistical  Methods,  Chapter  3,  McGraw  Hill,  N.Y. 
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where  D  =  oxxa yy  -  (oXy)~  and  0  ‘s  a  constant.  If  we  approximate  Ax  and  Ay  by  Gaussian 
random  variables,  then  the  probability  that  the  target  solution  lies  within  the  ellipse  (B-5)  is 
To  show  this,  we  rotate  the  coordinate  system  so  that  Ax  and  Ay  are  replaced  by 
variables  i  ]  and  i  ->  along  the  principal  axes  of  the  ellipse,  liquation  ( B-5 )  becomes 


-I  s 


O]  or; 


and  the  joint  probability  density  is 


[~2  2 

A  /  zi  z~> 

By  changing  to  polar  coordinates  and  integrating  (B-7)  with  r  /  — ±j-  +  -~ 

"V  °T  °2 

from  r=0  to  r=7ci  we  find  the  probability  to  be  l-e-0'-.  Typically,  c  is  taken  to  be  1  which 
yields  a  confidence  region  of  bV'r. 

COMPUTATION  OF  Q 

Equation  (8)  provides  an  expression  for  the  covariance  matrix  of  the  state  variables 
arising  from  the  measurements  of  a  single  random  variable.  For  several  such  independent 
variables  wj,  wi.  . .  .  wM,  the  linearization  (4)  becomes  a  sum  over  M  terms.  Since  these 
errors  are  independent,  the  covariance  matrix  resulting  from  their  simultaneous  presence  is 


°T-  I 


(B-8) 


where  0  (o“,. )  is  given  by  (8)  with  o^,  replaced  by  a~v,.  In  particular,  the  variance  of  the 

state  variables  is  (approximately)  linear  in  the  variances  of  the  measured  parameters.  The 
standard  deviation  is  the  root  mean  square  of  the  individual  contributions. 


A  plot  produced  by  the  error  estimation  software  appears  in  Figure  B-  la.  Although  not 
shown  in  that  figure,  the  error  covariance  QPP  is  a  decreasing  matrix  function  of  the  data 
(i.e.,  the  corresponding  error  ellipses  “decrease”  with  the  addition  of  each  data  point).  As 
previously  mentioned,  this  is  a  necessary  consequence  of  the  fact  that  the  least-squares 
weights  were  chosen  inversely  proportional  to  the  true  variances.  On  the  other  hand,  we 
note  that  the  ellipses  plotted  in  Figure  B-la  are  not  always  “decreasing.”  Such  behavior 
results  from  the  propagation  of  the  error  along  the  target  track;  the  covariance  at  a  given 
time  being  a  complex  function  (equation  (B-3))  of  the  velocity  covariances  QPV  and  Qvv 
as  well  as  QPP.  In  addition,  one  can  expect  a  discontinuity  between  the  two  target  segments 
since,  at  the  time  of  maneuver,  the  number  of  state  variables  increases  by  two. 
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Also  included  in  the  plot  are  points  computed  by  an  MLH  tracking  routine  from 
simulated  data.  The  weights  at  and  bt  of  equation  ( 1 ).  used  in  the  tracker  were  chosen  to  be 

(O0)~~  and  (of/f0)-“i  respectively.  Typically  these  are  a  priori  variances  and  hence  not 
necessarily  very  accurate.  However,  the  error  model  is  valid  for  any  choice  of  these  weights, 
independent  of  the  measurement  statistics;  those  used  in  the  model  should  simply  be  the 
same  as  those  in  the  cost  function  of  the  tracker.  On  the  other  hand,  the  predicted  errors  still 
depend  critically  on  the  statistics  through  the  variances  o^.  of  equation  ( B-8).  In  practical 
situations  these  must  be  estimated  from  the  data.  For  the  tracker  used  here,  the  variances 
for  the  bearings  and  frequencies  were  determined  by  averaging  the  errors,  and 

fj -f  j  with  and  fj  determined  by  the  solution  ( not  true)  track.  The  averaging  also 
took  into  account  the  reduced  number  of  degrees  of  freedom  resulting  from  the  dependence 
of  the  and  fj  forced  by  equation  ( 10).  The  results  were  used  to  plot  the  error  ellipses 

in  Figure  B-lb.  (Those  in  Figure  B-la  used  the  true  track  and  true  o^,.,  i.e..  og  and  oj: .) 

Note  that  the  o^,.  estimates  for  the  southern-most  ellipse  were  apprently  too  small.  This  is 

not  surprising  considering  the  small  number  of  data  points  (three)  involved  in  the  averaging. 
No  solution  was  found  for  less  than  three  data  points  since  the  number  of  state  variables 
(five  before  the  maneuver)  outnumbered  the  number  of  equations  (  2  x  number  of  data 
points  for  frequency  and  bearing). 

As  a  second  example,  Figures  B-2a  and  B-2b  illustrate  the  effects  of  an  uncertainty 
in  platform  velocity.  This  case  is  identical  to  Figures  B-la  and  B-lb  except  that  the  bearing  and 
frequency  measurement  errors  have  been  set  to  zero,  while  errors  in  platform  velocity  have 

been  introduced.  In  particular,  o'  =  var(v  )  (not  of  v'  !).  and  o“,  =  var(v,,)  with  the  values 

w  j  a  —  a  w  ~)  y 

chosen  to  be  =0"^  =0-5(knot)"  so  that  the  RMS  speed  error  is  1.0  knot.  These  errors 

correspond  to  errors  in  the  parameters  z  of  equation  (9)  rather  than  in  c  as  in  Figure  B-l . 

Also,  note  that  the  least  squares  weights  (taken  the  same  as  in  Figure  B-l )  are  not 
equal  to  the  inverses  of  and  /af/f0\  ",  both  of  which  are  zero  in  the  current  sample. 


Initial 

time 

(minutes) 

Final 

time 

(minutes) 

Course 

(degrees) 

Speed 

(knots) 

Target  leg  1 

1  ° 

33.33 

315.0 

|v,  1=15.0 

Target  leg  2 

33.33 

!  72.0 

10.0 

|v':  |  =  !5.0 

Platform  1 

0 

72.0 

30.0 

6.0 

Platform  2 

|  4.0 

1  64.0 

135.0 

5.0 

State  Variables:  P'QX,  P^,  v', x,  v,y,  v,x,  v,y,  fQ 


Table  for  Figures  B-)a  and  B-lb. 


Figure  B-la.  Comparison  of  predicted  errors  with  actual  results 
from  an  MLE  tracking  routine.  The  bearing  standard  deviation  is 
0/)=  1.5°  and  frequency  standard  deviation  is  af=  0.01  with  a  base 

v  V  /  » 

frequency  ot  100  Hz.  The  MLE  weightsare  a and  )”*"  - 
The  error  ellipses,  etc.,  were  generated  by  the  sensitivity  error 
analysis  routine  from  the  true  track.  (See  preceding  table  for 
added  information.) 


Figure  B-lb.  The  “reverse”  of  Figure  B-Ia.  This  plot  was  generated  using  an 
actual  MLE  tracking  routine  on  simulated  data.  Pseudo  random  noise  with 
standard  deviations  of  1.5°  and  0.01  Hz, were  added  to  the  true  bearing  and 
frequency  values.  The  error  analysis  producing  the  ellipses  used  estimates  of 
measurement  errors  and  target  positions  supplied  by  the  MLE  algorithm. 

(See  preceding  table  for  added  information.) 
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Figure  B-2b.  Same  as  B-2a  but  using  an  MLE  tracking  routine  with 
simulated  data. 


APPENDIX  C 

SOME  THEOREMS  ON  POSITIVE  DEFINITE  MATRICES 


THEOREM  I 

If  A  and  B  arc  positive  definite  symmetric  matrices  satisfying  A  >  B,  then  B“*  >  A. 
The  proof  relies  on  the  following  lemma: 

Lemma  I  ( |9| ,  pg  3071 

If  A  and  B  are  symmetric  matrices  with  A  positive  definite,  then  there  exists  an 

A  4- 

(more  than  one)  invertible  matrix  M  such  that  H  AH  and  H  BH  are  diagonal  and  real  (  + 
indicates  transpose). 

pf  of  Theorem  I 

Given  any  vector  x.  let  y  =  Hx.  Then  y '  Ay  >  y '  By  implies 

xfH+AHx  >  x+HfBHx  for  all  x  .  (C-l) 

Let  DA  and  Dg  be  the  diagonal  matrices  of  Lemma  I.  (C-l )  implies 

xTDAx>x'DBx  all  x  .  (C-2) 

Since  the  entries  of  DA  and  Dg  must  be  positive. 

x'D-a'x  <  xtDg1x 


or 

x+H~1A'1H+_,x<xi'H"1B"1H';  ~1x  for  all  x  .  (C-3) 

Given  any  z.  let  x  =  H+z,  then 

zf  A-1  z<  z+B_lz  all  z  ,  (C-4) 

which  proves  the  theorem. 

THEOREM  II 

Let  A  and  B  satisfy  the  hypotheses  of  Theorem  I  .  then  their  largest  and  smallest 
eigenvalues  satisfy  X^ax  >  Ag*ax  and  A™'11  >  Xg1  'n  . 


9.  Parlett,  B.,  The  Symmetric  Eigenvalue  Problem. Prentice  Hall.  In...,  N.J..  p  307.  1080. 
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pf  of  Theorem  II 


Xmax  = 

max 

llx||=l 

xAx  > 

max 

llxll 

xBx  =  \gax 

( C-5 ) 

xmin  = 

min 

llx||=l 

xAx  > 

min 

11x11=1 

xBx  =  Xg1*11 

(C-6) 

THEOREM  III 

Let  A  and  B  be  as  in  Theorem  I.  Then  for  any  given  c  >  0  the  concentration  ellipsoid 
Ea  =  ^x  |  x^  Ax  <  c|  lies  within  the  ellipsoid  Eg  =  jx  |  x+  Bx  <  c  j  . 

pf  of  Theorem  III 

Let  x  lie  on  the  surface  of  ellipsoid  E^.  Then  c'  =  x+  Bx  <  x^  Ax  =  c.  Let  r  =>/c/c' 

>  1.  Then  (rx)^B(rx)  =  c.  Thus  rx  lies  on  the  ellipsoid  Eg.  Then  by  convexity  x  e  Eg  and 
hence  E^  C  Eg. 
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